%--------------------------------------------------------------------------
function [x0, tf] = dro_(x0,auxdata)
mu = auxdata.mu;
RelTol = auxdata.RelTol;
AbsTol = auxdata.AbsTol;

MaxNum = 100; % 最大迭代次数
Maxtf = 10; % 最长积分时间
DAMP = 0.5; % 增加DAMP改善收敛性

% -------------------------------------------------------------------------
xf = ones(6,1);
count = 1;		% 开始计数

while (abs(xf(1) - x0(1)) > 1e-14 && abs(xf(4)) > 1e-14 && abs(xf(5) - x0(5)) > 1e-14)
    % 超出最大迭代次数
    if count > MaxNum
        disp('Maximum iterations exceeded') ;
        break
    end
    
    % 状态变量和STM
    xx0 = zeros(42,1);
    xx0(1:6) = x0;
    xx0(7:42) = reshape(eye(6) , 36 , 1);

    %求解ode
    options = odeset('RelTol',RelTol,'AbsTol',AbsTol,'Events',@event_);
    [tt,xx,tf,xxf,ie] = ode45(@fun_crtbp_stm_,[0,Maxtf],xx0,options,auxdata);
    xf = xxf(2,1:6); % 终止状态变量
    phif = reshape(xxf(2,7 : end),6,6); % 单值矩阵

    %更新相对距离
    d2 = (xf(1) + mu)^2 + xf(2)^2 + xf(3)^2;
    r2 = (xf(1) - (1 - mu))^2 + xf(2)^2 + xf(3)^2;
    d3 = d2^1.5;
    r3 = r2^1.5;
    
    %更新势能
    Ux = xf(1)  - (1 - mu) * (xf(1) + mu) / d3 - mu * (xf(1) - (1 - mu)) / r3;
    Uy = xf(2) - (1 - mu) * xf(2) / d3 - mu*xf(2) / r3;
    
    %dxf和yf的二阶导
    ddxf = Ux + 2*xf(5);
    ddyf = Uy - 2*xf(4);
    
    %求解初始修正量
    %delta_x0 = pinv([phif(1,1)-1, phif(1,5) , xf(4);
                     %phif(4,1) , phif(4,5) , ddxf;
                     %phif(5,1) , phif(5,5)-1 , ddyf])*[xf(1)-x0(1) ; xf(4) ; xf(5)-x0(5)];
    delta_x0 = pinv([phif(1,1)-1, phif(1,5);
                     phif(5,1) , phif(5,5)-1;])*[xf(1)-x0(1) ; xf(5)-x0(5)];
                    
    x0(1) = x0(1) - (1 - DAMP^count) * delta_x0(1);
    x0(5) = x0(5) - (1 - DAMP^count) * delta_x0(2);
    %tf = tf - delta_x0(3);
    
    count = count + 1 ;
end

dd = [];

end
% -------------------------------------------------------------------------

% 积分终止条件，即积分到x-z平面终止
function [value,isterminal,direction] = event_(t , xx , auxdata)
value = xx(2); % x-z plane
isterminal = 1;
direction = 1; % direction of event
end